Home > sgwt_toolbox > sgwt_cheby_op.m

sgwt_cheby_op

PURPOSE ^

sgwt_cheby_op : Chebyshev polynomial of Laplacian applied to vector

SYNOPSIS ^

function r=sgwt_cheby_op(f,L,c,arange)

DESCRIPTION ^

 sgwt_cheby_op : Chebyshev polynomial of Laplacian applied to vector

 function r=sgwt_cheby_op(f,L,c,arange)

 Compute (possibly multiple) polynomials of laplacian (in Chebyshev
 basis) applied to input.

 Coefficients for multiple polynomials may be passed as a cell array. This is
 equivalent to setting
 r{1}=sgwt_cheby_op(f,L,c{1},arange);
 r{2}=sgwt_cheby_op(f,L,c{2},arange);
 ...
 
 but is more efficient as the Chebyshev polynomials of L applied
 to f can be computed once and shared.

 Inputs:
 f- input vector
 L - graph laplacian (should be sparse)
 c - Chebyshev coefficients. If c is a plain array, then they are
     coefficients for a single polynomial. If c is a cell array, 
     then it contains coefficients for multiple polynomials, such 
     that c{j}(1+k) is k'th Chebyshev coefficient the j'th polynomial.
 arange - interval of approximation

 Outputs:
 r - result. If c is cell array, r will be cell array of vectors
     size of f. If c is a plain array, r will be a vector the size
     of f.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % sgwt_cheby_op : Chebyshev polynomial of Laplacian applied to vector
0002 %
0003 % function r=sgwt_cheby_op(f,L,c,arange)
0004 %
0005 % Compute (possibly multiple) polynomials of laplacian (in Chebyshev
0006 % basis) applied to input.
0007 %
0008 % Coefficients for multiple polynomials may be passed as a cell array. This is
0009 % equivalent to setting
0010 % r{1}=sgwt_cheby_op(f,L,c{1},arange);
0011 % r{2}=sgwt_cheby_op(f,L,c{2},arange);
0012 % ...
0013 %
0014 % but is more efficient as the Chebyshev polynomials of L applied
0015 % to f can be computed once and shared.
0016 %
0017 % Inputs:
0018 % f- input vector
0019 % L - graph laplacian (should be sparse)
0020 % c - Chebyshev coefficients. If c is a plain array, then they are
0021 %     coefficients for a single polynomial. If c is a cell array,
0022 %     then it contains coefficients for multiple polynomials, such
0023 %     that c{j}(1+k) is k'th Chebyshev coefficient the j'th polynomial.
0024 % arange - interval of approximation
0025 %
0026 % Outputs:
0027 % r - result. If c is cell array, r will be cell array of vectors
0028 %     size of f. If c is a plain array, r will be a vector the size
0029 %     of f.
0030 
0031 % This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
0032 % Copyright (C) 2010, David K. Hammond.
0033 %
0034 % The SGWT toolbox is free software: you can redistribute it and/or modify
0035 % it under the terms of the GNU General Public License as published by
0036 % the Free Software Foundation, either version 3 of the License, or
0037 % (at your option) any later version.
0038 %
0039 % The SGWT toolbox is distributed in the hope that it will be useful,
0040 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0041 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0042 % GNU General Public License for more details.
0043 %
0044 % You should have received a copy of the GNU General Public License
0045 % along with the SGWT toolbox.  If not, see <http://www.gnu.org/licenses/>.
0046 
0047 function r=sgwt_cheby_op(f,L,c,arange)
0048 
0049 if ~iscell(c)
0050   r=sgwt_cheby_op(f,L,{c},arange);
0051   r=r{1};
0052   return;
0053 end
0054 
0055 Nscales=numel(c);
0056 M=zeros(size(Nscales));
0057 for j=1:Nscales
0058     M(j)=numel(c{j});
0059 end
0060 assert(all(M>=2));
0061 
0062 maxM=max(M);
0063 %Twf_new = T_j(L) f
0064 %Twf_cur T_{j-1}(L) f
0065 %TWf_old T_{j-2}(L) f
0066 
0067 a1=(arange(2)-arange(1))/2;
0068 a2=(arange(2)+arange(1))/2;
0069 
0070 Twf_old=f; %j=0;
0071 Twf_cur=(L*f-a2*f)/a1; % j=1;
0072 for j=1:Nscales
0073     r{j}=.5*c{j}(1)*Twf_old + c{j}(2)*Twf_cur;
0074 end
0075 
0076 for k=2:maxM
0077     Twf_new = (2/a1)*(L*Twf_cur-a2*Twf_cur)-Twf_old;
0078     for j=1:Nscales
0079         if 1+k<=M(j)
0080             r{j}=r{j}+c{j}(k+1)*Twf_new;
0081         end
0082     end
0083     Twf_old=Twf_cur;
0084     Twf_cur=Twf_new;
0085 end

Generated on Wed 13-Oct-2010 13:36:39 by m2html © 2003